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1 INTRODUCTION 

This is a report for the NASA Consortium Program NCC2-51 16, entitled 
‘Methods for In-Flight Robustness Evaluation’. The goal of this program was to combine 
modem control concepts with new identification techniques to develop a comprehensive 
package for estimation of ‘robust flutter boundaries’ based on experimental data. The 
goal was to use flight data, combined with a fundamental physical understanding of 
flutter dynamics, to generate a prediction of flutter speed and an estimate of the accuracy 
of the prediction. 

This report is organized as follows: the specific contributions of this project will be listed 
first. Then, the problem under study will be stated and the general approach will be 
outlined. Third, the specific system under study (F-18 SRA) will be described and a 
preliminary data analysis will be performed. Then, the various steps of the flutter 
boundary determination will be outlined and applied to the F-18 SRA data and others. 


2 Specific Contributions 

The specific contributions of this project include development of a robustness problem 
formulation for flutter clearance, software development, and testing of general 
identification methodologies on recorded flight data. The software developments include 

• A new subspace identification algorithm that takes into account the presence of 
multiple data sets, 

This research was supported by NASA Ames University Consortium, 

Re: NCC2-5 116. 



• A novel algorithm to identify transfer functions in very noisy environments, 

• An algorithm to identify a parametric model of flutter dynamics based on frequency 
response using a Newton-based optimization procedure. 

• A procedure to project validated models to a predicted flutter boundary, based on 
modem robustness analysis methods. 

• A first-cut demonstration of the flutter boundary prediction procedure for NASA 
Langley wind tunnel flutter test data. 

In addition to the demonstration on Langley data, all of these algorithms have been 
tested on flight data provided by Dryden. Their validation is currently under study. The 
contributions related to subspace identification have all been reported in a paper that 
appeared at the AIAA Guidance, Navigation and Control Conference [1] and is provided 
in an Appendix. 

3 Problem Statement and approach 
3.1 Problem Statement 

Current flutter clearance procedures rely on post-processing of flight test data, 
which significantly slows the envelope expansion process for new or modified aircraft. A 
priori prediction of the flutter boundary is difficult, and flight test data are currently not 
used to systematically update these a priori predictions. Further, many methods for 
determining proximity to the flutter boundary are based on either da m ping ratio 
identification or tracking the evolution of spectral peaks. Experience and analysis have 
shown that these measures can be highly nonlinear with respect to flight parameters such 
as dynamic pressure and Mach number, and thus may be misleading measures of flutter 
stability (we present a simple illustration of this problem below). A method is needed to 
alleviate these difficulties, that systematically combines a priori knowledge on the 
aircraft with flight test data into a well-behaved prediction of the flutter boundary, and 
that does so in near real time so that envelope expansion can proceed more efficiently. 

Current practice of using damping ratio as a measure of proximity to instability is 
motivated by the fact that as an aeroelastic mode approaches neutral stability, its 
damping ratio approaches zero. Thus one way to predict the flutter boundary is by 
tracking the damping ratio of each flexible mode as a function of critical parameters 
(dynamic pressure or Mach number), and subsequently predicting the flutter boundary by 
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simple linear extrapolation. In many cases however, such a measure can be unreliable. 
For example, we have plotted in Fig. 1 the evolution of the damping ratio as a function of 
speed for a typical wing section [8]. We see that a linear extrapolation of the damping 
ratio will systematically overestimate the flutter boundary. Besides, when far away from 
the flutter boundary, the damping ratio actually increases as the aircraft approaches 
flutter. 


TypiolStctioa Problem 



Damping ratio versus speed 



Figure 1 - Evolution of damping ratio vs. speed for a typical section. 


32 Approach 

The approach proposed in this research may be sumarized in the following 
diagram: 



Subspace ID or Nonlinear Robustness 

direct TF estimation least-squares Analysis 


To develop a reliable flutter clearance capability, we cast the flutter clearance 
problem as a ‘robustness problem* typical in control applications. The flutter clearance 
problem stated as a robustness problem is simply this: based on the currently available 
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data and models, what is the ‘smallest’ perturbation to the current flight condition which 
will drive it to the flutter boundary? Size of the flight condition perturbation is 
represented as a norm or vectoral distance; this type of measure is typically much more 
well-behaved than damping ratio or spectral peak amplitude. It is called the structured- 
singular value in control theory, and it is the multivariable counterpart of the classical 
gain margin for single-input, single-output control systems. Several tools are now 
available that compute it quickly and accurately. 

To incorporate flight test data into the robustness problem, we take a unique 
approach which involves compressing information from multiple data sets, taken at 
various flight conditions. This is necessitated by the quantity and quality of the data 
available (large amounts of relatively noisy data), as well as the nonlinearity of the 
variation of the dynamics with flight condition. Model structure information is retained 
by the procedure, through a nonlinear least-squares fit of the system dynamics to the 
acquired data. 

To ‘compactify’ the data and allow more efficient identification, subspace 
identification methods have been investigated (see for instance De Moor and Vandewalle 
[1] and Van Overschee and De Moor [2]). Alternatively, we also have investigated the 
possibility of compressing data by directly generating reliable, de-noised transfer 
functions. Both technique are essentially parameter-free identification procedures, which 
provide models that are much reduced in size compared to the original time domain data, 
but which are completely unstructured and therefore require additional identification to 
be useful for robustness analysis. This additional identification is performed by adjusting 
optimally the coefficients in an aeroelastic model that is known a priori, such that both 
models agree as well as possible from an input-output viewpoint Currently the 
optimization criterion is a frequency weighted H2 norm. 

4 System description and preliminary data analysis 

4.1 System description 

The F-18 SRA (Systems Research Aircraft) is a test aircraft that has been equipped with 
sophisticated flutter envelope clearance equipment, including a pair of aerodynamic 
wingtip exciters that provide oscillatory wingtip lift. Ten accelerometers are located on 
all the aircraft’s body to measure the aeroelastic effects of the excitations. The table 
below summarizes the numbering and position of each sensor. 
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Output 

sensor location 

i 

left wing forwarc 

2 

“ left wing aft 

3 

left aileron 

4 

left vertical tail 

5 

left horizontal tai 

6 

right vertical tail 

7 

right horizontal te 

8 

right wing forwar 

9 

right wing aft 

10 

right aileron 


Accelerometer locations 


The inputs were chosen to be sinusoidal sweeps spanning the 3-30Hz range. This range 
was chosen because it is expected to contain all the flutter modes for the F-18. The 
sampling frequency was chosen to be 200Hz. Each sweep lasted 30 seconds to 
compromise between the need for reliable information and the requirement to save on 
fuel and maintenance costs. 

To be sure to excite the symmetric and antisymmetric dynamics of the aircraft, each test 
consisted of two sequences. In the first sequence, the two exciters were roughly in phase, 
whereas in the second experiment, these exciters were roughly 180 degrees out of phase. 
The tests were performed at different elevations (10K, 30K and 40K feet) and different 
Mach numbers (.8, .85, .9, .95) to obtain a broad range of flight conditions. The 
operating conditions at which tests were performed is plotted on the graph below, along 
with the aircraft flight envelope and the assumed flutter boundary. 
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42 Data analysis 


The data were collected from a real flight test and were significantly corrupted by 
atmospheric and other disturbances. In addition, some sensors were suspected to be 
defective. To investigate both issues, coherence plots between each input and each output 
were computed. The corresponding coherence plots may are illustrated in the figure 
below, which is based on a flight test performed at Mach 0.8 and 10K feet. 

It may be immediately remarked that on average, the measured coherence is low 
(no more than 0.8 in most cases). This indicates that the data at hand are contaminated by 
high levels of noise. From experiment to experiment, the coherence was also found to 
change significantly (possibly due to different weather conditions); such a difference in 
coherence may be used to weight results from many experiments differently. To quantify 
the results that have been obtained, a norm on the coherence plots was chosen to be the 
mean of the coherence between 0 and 50Hz, bearing in mind that the excitation 
frequency range is from 3 to 30 Hz. The resulting norms are tabulated below. 



Co fl£ &ev<-c 
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Coherence with left input 
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Flight condition for each data set 


The average of this norm for each outputs and for both inputs was calculated; it appears 
that five outputs consistently have a much higher score than any other outputs. These five 
outputs are the leading and trailing edge accelerators on each wing, and the right aileron 
accelerator. It was decided to discard the latter, because it also measures the aileron’s 
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own dynamics, which are neglected in the ensuing analysis. When more experience is 
acquired, this viewpoint may be revised. 

Browsing through the public domain literature has led to interesting comparisons. For 
example, Bucharles, Cassa and Robertier [7] use signals with average coherences ranging 
from 0.9 to 0.95 to perform the flutter analysis of Airbus commercial aircraft. 

5 Data compression 

Two data compression approaches have been investigated during this research. First, 
subspace identification algorithms have been tested and improved. Second, because the 
data under consideration are so noisy, an alternative direct transfer function evaluation 
scheme was developed, which is robust to a large range of noise disturbances. 

5.1 Subspace identification: review of existing techniques and new developments. 

As most of the information gathered during this phase of the project has been already 
reported in the conference paper [1], this section will only highlight its main findings. 

The attached appendix contains the conference paper. 

Subspace identification techniques represent a way to perform system identification in an 
unstructured way, using state-of-the-art results in matrix computations, including least- 
squares and total least-squares techniques. Subspace identification techniques are 
essentially “black-box” algorithms, which attempt to infer a low-order, state-space 
model that “best” explains a given set of inputs and outputs. To evaluate how useful 
these procedures are for aerospace applications such as flutter boundary determination, 
the F-18 SRA data were used as benchmark data. Many subspace algorithms currently 
available were tested on these data, especially the algorithms by Cho and also by 
DeMoor [l]and VanOverschee [2]. 

The results of this study led to the following conclusions: 


(i) The collected data are often distributed over many flights at the same flight condition 
This is incompatible with existing subspace identification software, which were all built 
for single I/O data sequences. To remedy this problem, the basic principles of subspace 
identification were reviewed and adapted to handle many, uncorrelated data sets. This 
technique has since then found applications in other areas of interest to NAS A-Diyden, 
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including the reliable identification of the F- 18 SRA lateral-directional dynamics for 
nonlinear simulation calibration and autopilot design purposes. 

(ii) The data are very noisy and subspace identification methods are on the whole unable 
to cope with these levels of noise. Besides, different subspace identification techniques 
do not all share the same performance. As a result, an alternate data compression scheme 
whereby transfer functions are directly estimated has been devised. 

52 Direct transfer function estimation 

This method is based on the a priori knowledge that the input signal is a frequency 
sweep. Thus it is possible to predict a priori at what time the relevant frequencies should 
appear in the input and output signals. 

5.2.1 Single input-single output case 

We first consider the single input-single output case. The assumptions for this method is 
that the excitation is a frequency sweep. Process and observation noises are present as 
shown in the following figure. 


y 

Input 

noise 

k J 

System 

Gif co) 

^ 

Output 

noise 

V ' 

Input 

signal 

rv - >/ Output 
signal 


The step by step procedure of the transfer function estimation is described as follow. 

Step 1: localization of the relevant frequencies by filtering the signal with a set of narrow 
band pass filter centered around the frequencies C0/ . Filtering the input and the output 
does not alter the identification. Indeed if the input and the output of the system are 
filtered with the same linear time invariant filter F, the filtered output appears as the 
output of the system excited by the filtered input. 
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Step 2: localization of the relevant time information by windowing each filtered signals. 
The chosen window has to be the same for the input and the output and is chosen to be 
centered at the point where the filtered input has the maximum amplitude. At this step, 
the two signals should appear as sine functions. 


Step 3: calculation of the estimate at the frequency by the ratio of the output Fourier 
transform of the windowed signal to the input one: 
y V ln)e 2 * j< °' n 

G( /co ) = r v' 7ie * - — — — , where -^ mes and ^“denotes the measured signals. 

Let us now look at the properties of this estimate. Let ^ and ^ be the input and output 

noises, respectively. Then we have 

, G(yco,.)(Viyco / ) + /^(M)) + A4(yco ( .) 

G(yco ' } m W) 

-re \ l rtt \ , ^(M) 

-G(yco,) + G(yco ( ) 


N ( i(£> ) = y 77 ( k)^ k = 0 

If the two noises are assumed to be unbiased, then ' w '' ,v ' 

Therefore, 6( yico ( ) = G(yco, ) , i.e. the estimate is unbiased. 

The variance of the noise is given by G N ~ ) - N 0 a > w here ° is the 

variance of the noise, and is the number of points in the time window. The variance 
of the transfer function estimate can now be calculated, assuming that the two noises are 
uncorrelated: 


<*g 2 = 


G(yco ( ) 


U{ M) 


N 0 o* + 
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U( M 




No 


UUtoi)=?U 

Assuming that the windowed input is a sine function, 2 , with U being the 

amplitude of the input, the variance of the estimate is 

, . ,2 4a. 2 4a, 2 

I ■ 


N 0 U 2 


NqU 


5.2.2 Multiple input single output case 

At this point, let us assumed that the system has m inputs. A set of m independent 
experiments is necessary to use the proposed method and all the excitations are assumed 
to be frequency sweeps. The first step is identical to the single input case. The property 
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that the filtered output will be identical to the output associated with the filtered inputs is 
still satisfied since, in the Fourier domain, the matrix of the filter is a scalar matrix which 
commutes with any other matrix. For the second step, the procedure is the same but one 
of the m inputs has to be taken as a reference in order to select the time window. For the 
third step, the Fourier transform of every signal has to be calculated and the estimate is 
given by 


G|(M) 


X(M) 


-1 







YJL M)_ 


with ^ k ' the 1 th input of the k* experiment. Note that the input matrix has to be full rank 
in order to apply this formula (if it is not full rank, we conclude that t here not enough 
information in the data to estimate the transfer function). With the same assumption for 
the noise, the estimate still has no bias. The variance of the estimate is given by 

f 3(M)-§(M)1 



[<$ (M)~ Gt ( M ) ••• (j®i) -G\(j®m) } 




where E stands for the expected value. By plugging in the value of the estimate, this 
becomes 
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Since all the noises are assumed independent, this equation becomes 
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Since the middle matrix is a scalar matrix (diagonal with all its diagonal terms equal), it 
commutes with the input matrix. Thus the variance is proportional to the singular values 
of the inverse of the input matrix. This matrix is independent of the noise. Its singular 
values give some information on how well the inputs have been chosen for identification 
of the system. 


5.2.3 Numerical experiment 


This method was applied to simulation data using a model of the structural dynamics of 
the FI 8. The system has two inputs. The two sets of experiments needed for the 
estimation were symmetric and asymmetric excitation, with equal amplitude U for each 
input. For each time k, the process noise n(k) was chosen to be a random variable 
uniformly distributed over the interval [-U/4 U!4\, thus the variance of each of the inputs 
is U 2 /(4*12). No sensor noise was added in this experiment. The time window length 
was chosen to be 128 points. The variance was estimated with the simulated data by 

using the classical estimator for the variance, ^ ** , where N is the number of 

experiments and x is the error of the estimate. One hundred simulations were made in 
order to obtain a 90% confidence interval for the variance of approximately 5%. The 


input matrix can be written as follow: 
theoretical variance of the estimate is 
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J. Since the output noise is 0, the 
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Note the variance of the transfer function estimate is independent of the noise amplitude. 
The variance was compared with the theoretical results and plotted below. In this case, 
the standard deviation of the variance estimate is 5% of the estimate and is independent 


of the frequencv. 

iruutftr (unci ton viriiinc. oi lh« trantler lunctwn 



5.2.4 Application to F18-SRA data 

This method was then applied on the real flight data using only the four outputs chosen 
in Section 2. The four transfer functions plotted below were obtained from F-18 SRA 
flight data taken at 10,000 feet and Mach 0.8. The first transfer function (input 1 to 
output 1) starts with a slope of 40 dB/decade at low frequency and a phase of 180 
degrees. The phase then drops to -360 degrees and the magnitude is constant, at around 
40 dB. This means that there is a second order pole at 6.5 Hz. It seems that there is also a 
pole zero cancellation around 13 Hz. However, due to the low resolution of the estimate, 
this cannot be affirmed by this plot only. The high frequencies (above 18 Hz) are rather 
noisy, but it seems that the magnitude and the phase are stable, meaning no pole or zeros 
are in this region. Looking at the second transfer function (input 1 to output 2), it seems 
that the low frequencies have the same properties as the first plot which is a slope of +40 
dB/decade for the magnitude with a phase of +180 (which is the same as -180). The pole 
at 6.5 Hz is still detectable on this plot and a pole around 12 Hz is also seen. This 
confirms the pole zero cancellation of the previous plot. For the asymmetric excitation, 
two poles could be detected by looking at the magnitude of the transfer function at 8 and 
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18 Hz. However, it is much harder to correlate the phase in this case, because its 
variation with frequency is much smoother. 



6 IDENTIFICATION OF STATE-SPACE DYNAMICS BASED ON COMPRESSED DATA 

Once the data are compactified (either in state-space via subspace identification or via 
direct transfer function estimation), it is necessary to fit a physical model of the system to 
them, 

where the input/output form of this model is given by 

\x= AX+ BU 
|V= CX+ DU 

The par ame ters in the fitting procedure are the significant elements of the system 
matrices (A, B, C, and D), where the system states are the physical states, and the initial 
guess and structure of the model (structure is represented by the zero and identity 
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elements of the system matrices) come from a priori modeling. In its basic version, the 
fitting operation consists of minimizing the H2 norm of the difference between the 
system defined by the matrices (A, B, C, and D) over the relevant coefficients (such as 
stability derivatives and aerodynamic lags) and the experimental, compressed data. In a 
more sophisticated version, a series of fits is done first over a range of flight conditions. 
These fits are then inteipolated using a specific parameterization with respect to, for 
example, Mach number and dynamic pressure q. 


6.1 Nonlinear fit via Newton algorithm 

To find a system of the form (1) that best matches the experimental data, we propose to 
constraint some elements of the matrices to be constant in accordance to the nominal 
flutter model that was chosen. The chosen cost function is the H2 norm of the difference 
between the estimated transfer function and the transfer function of the identified model. 
Even though some elements of the matrices were fixed, the optimization problem is still 
unconstrained in the mathematical sense since all the coefficients that are allowed to vary 
are totally free. 

The H2 norm for a system is given by: 

J= \Tr{(C(ja>- A)~'B+ D-G(yco))(C(yco- A)~'B+ D- G(yco))*)do 
0 

The gradient of J with respect to all the matrices needs to be evaluated. With respect to 
the matrix A, the gradient is: 

|^=2Re(JYr((C(y'co- A)~'dA(ja- A)~'B)(C{j(»- A)~' B+ D-G{jG>))')dto) 
dA o 

By using the multiplication commutation inside the trace operator, we have 

|^= 2Re(J 7r(((yco - A)~' B)(C(yco- A)~' B+D- G(yco))‘ C(yco - A)-'dA)dn) 
dA J Q 

With the same procedure, the gradient with respect to B, C and D can be evaluated: 

— = 2 Re(J Tr{ C{ yco - A)' 1 B+ D-G{ja))' C{ yco - A)~'dB) do) , 

dB J 0 

-^ = 2Re(f Tr((yco - A)" 1 B{C(j(i)- A)~' B+ D- G( ja))' dC)dco) , and 
dC J Q 

~\J 

— = 2 Re( J Tr{{ C{ yco - A) -1 B+D-G( yco ))‘dD) do) 
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These integrals are well approximated by a discretization (in frequency) at all the points 
at which the transfer function has been estimated. The gradient of the function J can then 
be derived from the above formulas. 

For a Newton optimization, it is also necessary to calculated the Hessian H of the 
function. However, this operation is computationally very intensive; it is generally more 
efficient to use a so-called quasi-Newton method instead. This method estimates H based 
on the variation of the value of the function and of the gradient. A efficient method to 
estimate H 1S the BFGS method. Once H is estimated, the direction of search is 
calculated by d=H~'VJ An optimization of J along this search direction is then 
performed to find the minimum J in the current search direction (i.e. the step length is 
optimized). The classical ID Newton method can be used at this stage, but it is then 
necessary to calculate the derivative in the direction of search at each step. Studies on 
optimization have shown that it is not necessary to find a very accurate minimum in the 
current search direction because, at the next iteration, the direction of search is going to 
change. The little gain obtained by optimizing in each search direction can usually be 
avoided if it is time consuming. The procedure that was used was developed by de Wolf. 
The derivative in the current search direction is estimated only at the starting point and 
noted P. Then a stopping region is defined by two coefficients ^ and ^ (these can be 
tuned but can be set to .4 and .6 respectively as a first try). ^ and ^ define two straight 
lines whose slope are respectively P^and P 17 ^. The region between the two lines is the 
‘stopping region’. A binary search is then performed until a point of the surface J falls 
inside the stopping region. 



The whole process is then iterated until a convergence is reached. The main drawback of 
this method is that it guaranties convergence to a local minimum only, and there is 
absolutely no way to distinguish, even a posteriori, if the point that was reached is a 
global or local optimum. 
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62 Application to F-18 SRA data 

6.2.1 Physical flutter model 

The first step to use the quasi Newton algorithm is to determine which elements in the 
state space model are fixed and which are to be optimized. Even though numerous flutter 
models have been proposed, a literature search showed that the most common is written 
as follows: 

0 / 0 

A= -M~'K - M-'C - M~'P 

0 0 A 

U 7 

where I is the identity matrix. The matrices M, K and C are the apparent structural mass 
stiffness and damping of the aircraft which means that that the non-circulatory part of the 
aerodynamics is included in the matrices. The lower-right matrix A represents the 
aerodynamic lags and is a diagonal matrix with real negative eigenvalues. P is the 
coupling term between the lags and the structure. 

6.2.2 Actual data and high-order model 

A finite element study of the structural dynamics of the FI 8 was computed by NASA 
Dryden. The full model was composed of 14 symmetric modes and 14 anti symmetric 
modes. 13 of them were in the 3-20Hz range and 11 of them may be involved in the 
flutter mechanism. It was therefore decided to reduce the model to 11 modes (or 22 
states). The first experiment that was tried was to put all the uncertainty into the lags. 
The matrices M, K and C were set by reducing the theoretical model to keep only the 22 
states of interest. For this, the state space system was diagonalized and only the modes 
with the appropriate eigenvalues were kept reduced. The matrices P and A could vary in 
the optimization and D was set to 0. It was assumed that the actuator acted as a force 
input on the structure and did not affect the aerodynamics around the wing. Therefore, 

T 0 " 

B= B 

the structure of the matrix B is chosen to be <-°J. The sensors are assumed to 
measure only structural displacements and no aerodynamics at all. In the model this 

shows in a matrix C of the form: ^ °1 Finally, it was found to be necessary 

to add an additional constraint in the Newton algorithm that forced the lag coefficients to 
stay above a certain limit (so that the identified model remains stable). 
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The optimization set up as described above saturates the constraints on the lags, and does 
not converge to even a local minimum. Thus it was decided to optimize over a single 
output, 2 outputs, and 3 outputs, to determine the effect of trying to fit more and more 
outputs on the optimization process. First, only one output was used and the algorithm 
ran with exactly the same structure of matrices. A stable local minimum which did not 
saturate the constraints was found for this case. A measure of the accuracy of the 
identification was defined to be the ratio between the value of the cost function J and the 
nominal H2 norm of the estimated transfer function, evaluated in %. This measure went 
down to 13.3% for the single-output test. When a second output was added, the accuracy 
of the identified model degraded to 20.6%. When a third output was added, the accuracy 
of the identified model stalled at 31.1%. Thus we conclude that it is progressively more 
difficult to explain all of the data with the single model that we are using; when all four 
outputs are chosen the search fails altogether. 

At this point, the order of the system to be identified was further reduced to 6 modes . 
Two things motivated this approach. First, computation time grows exponentially with 
both the number of parameters to fit and the order of the system, since it is necessary to 
calculate the inverse of a complex matrix of the size of A at every iteration. The second 
motivation was to detect whether some of the modes in the models were absolutely 
necessary. In this test, the system was modeled with only 6 modes. This approach 
converged with all four outputs, and yielded a system of order 18. The solution found 
had an accuracy of 33%. The figure on the following page shows the frequency domain 
plots of all 8 transfer functions, together with the results of the optimization. Despite the 
large value of the optimality measure, the basic shapes of the transfer functions have 
been captured consistently across all eight transfer functions. 

Based on these tests, we conclude that by reducing the model order, the difficulties 
encountered with convergence were alleviated to a sufficient degree that the four-output 
test was able to converge. The trend of error in the fit continues as expected, however 
the percent error is highest for the four-output identification. It is important to note, 
however, that the reduction in the model order did not significantly effect the overall 
accuracy of the fit (31% vs. 33%). Thus the hypothesis that some of the structural modes 
are indeed unnecessary to explain the data appears to be a valid one. Follow-on work is 
focusing on judicious model order reduction, as well as optimization of the mode shape 
parameters (contained in the C matrix) separately from the lag parameters. 
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6.2.3 Simulated data and low-order model 


Based on the results of the previous section, a test was executed using simulated data and 
a very low order system. As before this was motivated by the desire to reduce 
computation; for a very low order representation the number of local optimum should 
also be reduced. Therefore, since the Newton algorithm converges to any local minimum, 
the chances that the global minimum is reached are increasing. At this point, the physical 
model that was used primarily was totally ignored. The number of modes was set to two 
because two resonance appeared clearly on the transfer function estimates. The initial 

guess for the A matrix was of the following form: 

0 0 1 O' 

0 0 0 1 

A= , 


L 0 —l<2 C^2 ~ 

k t and k 2 were set so that their square root is equal to the approximate value that was read 
on the transfer function plot. The values of c were set so that the system is stable. These 
values were also chosen rather small since it is well known that structural dynamics 
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eigenvalues are lightly damped. The number of outputs was reduced to two because the 
sy mm etric and anti- symmetric part were assumed to be completely decoupled. The 
convergence of the Newton algorithm was quite fast as expected. More importantly, the 
fits were acceptable (less than 10% of relative error, see comparisons below), even 
though the model used for simulation had 14 modes and the model used for identification 
h ad only two modes. The experiment was tested for three different altitudes (10K, 30K 
and 40K feet) and the Mach number was kept constant at .8. At every flight point, both 
the symmetric and anti symmetric cases were tested. The results of the identification 
were always good and the convergence fast if a reasonable initial guess was utilized. 


irpt— >oul1 irp 1 — > out 2 




At this point, even if we assume that the optimum that was found was the global 
minimum, there is no guaranty that the three models are in the same basis. Therefore, 
some additional constraints were added to the second and third identification, forcing 
them to have the same C matrix as in the first identification. 

Another experiment was tried with the same set of data and the same parameterized 
model which consist in identifying the three Taylor expansion matrices in one single 
Newton algorithm. Since only three flight points were considered, the Taylor expansion 
should exactly match the previous results at the three elevations that were considered. 
The initial guess for this experiment was taken as the result of the previous Newton 
optimization for the nominal point, and the two other matrices of the Taylor expansion 
were set to 0. This means that the initial guess assume a constant A matrix with respect to 
the parameter. 
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6.2.4 Conclusions about F-18 SRA model identification 

Several critical issues were identified for the success of applying the Newton algorithm 
to identification of the structural dynamics of the FI 8. The computation time involved to 
solve the problem can be quite high if the model order and number of free parameters is 
not carefully considered. Choosing a reliable model order is a critical step, and further 
work is needed to determine what constraints should be applied and how knowledge of 
the physical system can be incorporated. Finally, the coefficients of the matrices should 
not be completely independent; at the very least they should be tied to dynamic pressure 
in a physically meaningful way - this is the focus of current research. 

6.3 Analysis of the Benchmark Active Control Technology (BACT) Demonstrator at 
NASA Langley 

Initial development of procedures for flutter boundary prediction were also tested 
in cooperation with NASA Langley, utilizing data from the BACT rig. This rig is a 32” 
wing section mounted in Langley’s transonic wind tunnel facility [3]. The purpose of the 
device is to test flutter control strategies in a benchmark environment. Because of this 
role, high frequency forcing can be introduced via trailing edge flaps and upper and 
lower surface spoilers. Spectral transfer function estimates of the aeroelastic behavior 
have been obtained, and detailed modeling has been conducted. Again, a nominal, finite 
element model of this system was available to be fit with these experimental, compact 
data, and we used the same Newton procedure to match this model with these data. In 
this case, a series of fits was done first over a range of flight conditions. These fits were 
then interpolated using a specific parameterization with respect to, for example, Mach 
number and q. The parameterization in use in this case is a quadratic of the form 

A = A q 2q 2 + A q iq + A^M 2 + A^iM + AMqMq + Aq. (*) 
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Experimental (curved line) and predicted (x’s) flutter boundaries for BACT using 
approximate aeroelastic model and experimental frequency-response data. Circles in 
plot at right correspond to flutter points (x’s) on plot at left. 


Typically, the matrices in the summation have the same structure as the original 
aeroelastic model. However, the parameterization in Mach number and q does not 
necessarily reflect the one originally present in the aeroelastic model. Rather, it is to be 
seen as a low-order approximation of it, which is flexible enough to incorporate badly 
modeled parameter dependencies such as Mach number. 

The plots above show the results obtained by applying this procedure to the 
Langley BACT data. The results are shown both on an M-q plane, and in terms of 
percent prediction accuracy. The latter plot requires some explanation. If the procedure 
is used during flight test, its accuracy becomes more important as the flight condition 
approaches the flutter boundary. Because the prediction is projecting to a boundary 
which is not as far away, the prediction also becomes more accurate as the flight 
condition approaches the flutter boundary. At any flight condition, as a minimum, the 
algorithm must say that it is safe to fly 50% of the remaining distance to the boundary. If 
it does not, then the envelope expansion will be halted prematurely. Overprediction of 
the flutter boundary may be acceptable when one is far from that boundary, but as it is 
approached the prediction must become more accurate. 

These considerations lead to the plot at left - if the predicted flutter speed falls 
within the shaded region at all times, then the prediction is considered to meet the 
minimum requirement for usefulness during envelope expansion. For the BACT it is 
seen that we are very nearly successful based on only 6 data sets (this is the minimum 
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number of data sets necessary to eventually find a parameterization of the form (*)). It 
remains to be seen if this accuracy can be obtained for data from a real aircraft 
Furthermore, since in the real case the flutter boundary will not be known, some measure 
of the accuracy of the prediction must be developed. 

7 REMAINING RESEARCH ISSUES 

Our attempts at developing the methodologies outlined above have raised the 
following issues: 

• As compared to the BACT data used at Langley, the F-18 SRA data is more 
complete; however it is corrupted by higher noise levels. As a consequence, 
existing subspace identification procedures have difficulty providing accurate, 
compact representations of the data. Improving identification accuracy will 
require research on modifications to subspace identification methods (for 
instance, frequency-dependent weighting). The relationship between the 
efficiency of the subspace identification procedure and the excitation signal at the 
wing tips is also an open research question. Discussions with members of the 
Structures group at Dryden suggest that some alternate excitation schemes might 
be implementable on existing hardware (currently logarithmic and linear sine 
sweeps are used). 

• The current nonlinear least-squares procedure to match the aeroelastic model with 
the ‘compactified’ data is computationally greedy, although specialized time- 
saving concepts such as Quasi-Newton have been used. While the size of the 
Langley system allowed the procedure to converge acceptably fast, improvements 
still need to be made for the same procedure to work on F-18 SRA aeroelastic 
models. 

• The choice of an appropriate parameterization as a function of Mach number and 
q is still not clearly defined. The parameterization (*) is such that after optimal 
fitting, some entries in the matrices A q 2 , A q i, Am 2 . Ami, AMq or A 0 may not be 
physical (in other words, the aeroelastic model, if used alone, would predict such 
entries to be zero). Schemes that impose the structure of the aeroelastic are 
possible; however, the 'unstructured' approach allows one to correct for badly 
modeled parameters such as Mach number. 

• In order to convert the flutter boundary problem into a robustness problem, 
various levels of sophistication are possible: in the most basic version, only the 
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Mach and q parameters appear in the perturbation block A. However, one can 
imagine that in the future, the two aforementioned steps (data compactification 
and least-squares model matching) provide parameterized models with error 
bounds; these could be easily and effectively inserted in the robustness problem 
as true additional uncertainties, leading to a more reliable analysis without 
necessarily introducing unwarranted conservatism. 

• In order to make the tools usable during flight test, a significant decrease in the 
computation time is required. Once methods with desirable properties have been 
validated, optimization of the computations must be considered. Recursive 
implementation, for instance, should be investigated. 

• An ongoing concern is proper interfacing of the proposed tool with the flight test 
engineer. In particular, we will need to carefully examine how data pertaining to 
flutter proximity should be displayed and documented in a reliable fashion. Also, 
any procedure able to detect malfunction of the proposed method (for example, 
bad fit of the experimental data) should be implemented and displayed as a 
warning to the operator. 
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